* This is the master do-file that runs all of the code necessary to replicate
* "Railroads, Reallocation, and the Rise of American Manufacturing" by
* Richard Hornbeck and Martin Rotemberg.

set more off
cap ssc install ereplace
cap ado uninstall reghdfe
net install reghdfe, from("https://raw.githubusercontent.com/sergiocorreia/reghdfe/master/src/")
cap ssc install carryforward
cap ssc install reclink
cap ssc install unique
cap ssc install moreobs
cap ado uninstall ftools
net install ftools, from("https://raw.githubusercontent.com/sergiocorreia/ftools/master/src/")
cap ssc install twostepweakiv
cap ssc install avar
cap ssc install ivreg2 
cap ssc install moremata 
cap ssc install ranktest
cap ado uninstall ivreghdfe
net install ivreghdfe, from("https://raw.githubusercontent.com/sergiocorreia/ivreghdfe/master/src/")
cap ssc install hdfe
ssc inst _gwtmean
cap log close

global flg_matlab = 1 // set to 1 to recalculate matlab parameters

*cd to proper directory
cd

if "`c(username)'" == "mr3019" {
	global main_path = "C:\Users\mr3019\Dropbox\MFG_Productivity\Final Data and Analysis"
}


if "`c(username)'" == "wcockrie" {
	global main_path = "C:/Users/wcockrie/Dropbox (UChicago)/MFG_Productivity/Final Data and Analysis"
} 

if "`c(username)'" == "abers" & "`c(os)'" == "Windows" {
	global main_path = "C:/Users/abers/Dropbox (UChicago)/RR_Replication/Final Data and Analysis"
} 

if "`c(username)'" == "abers" & "`c(os)'" == "Unix" {
	global main_path = "/project/cmf/rail/Final Data and Analysis"
} 

* This file uses Stat/Transfer to convert a dbf file to stata dta
* this program can be downloaded at: https://stattransfer.com/downloads/windows_downloads/
* This code requires Stat/Transfer 14 installed.
* set the value of StatTransfer_path to where st.exe lives. This is needed for command stcmd to work
global StatTransfer_path "C:/Program Files/StatTransfer16-64" 
* Define Conley Standard Error Program (note that ado file is included, but this do-file will accomplish
* the same goal on any system with any file structure
cd "$main_path"
cap program drop ols_spatial_HAC_miles
do "generate_prestructural/ols_spatial_HAC_miles.do"

********************************************************************************
***** GENERATE: PRE-STRUCTURAL *************************************************
********************************************************************************

cd "${main_path}/generate_prestructural"
do generate_prestructural.do

********************************************************************************
***** STRUCTURAL: OPTIMIZE THETA, COUNTERFACTUAL CALCULATTIONS *****************
********************************************************************************
cap program drop check_matlab
program define check_matlab
	glo file_is_missing=1
	while $file_is_missing {
		cap confirm file "intermediate/MATLAB_finished.csv"
		if _rc {
		* Stata checks if the output is complete approximately every second
			sleep 1000
		}
		else {
			erase "intermediate/MATLAB_finished.csv"
			glo file_is_missing=0
		}
	}
end

cd "${main_path}/structural"

if $flg_matlab==1 {

clear
	* Set calc_theta_full equal to 1 if you are reoptimizing theta
	local calc_theta_full = 1
	* Sets the wedges and cost shares to be used for theta/price optimization
	global robustness_str = "_data1880"

if `calc_theta_full' == 1 {

	* This path could take weeks depending on your computation power available
	* It will create several GBs of files, only one of which is truly needed
	* We define the range of thetas and avg prices that we want to test. The size
	* of each of these ranges will impact the time it takes to run dramatically.
	* The could be run in parallel if you are equipped to handle that with a
	* little bit of coding adjustment.
	global firsttheta = 1.5
	global lasttheta = 6.5
	global iters = 100
	global first_avgprice = 20
	global last_avgprice = 60
	global step_size = 1

	shell matlab -r "cd '${main_path}/structural';fst=${firsttheta};lst=${lasttheta};iters=${iters};avgprice_list=[${first_avgprice}:${step_size}:${last_avgprice}];run Data_and_Theta_Prep_gridsearch.m;exit"
	* This waits for the last file to be created, should be very long
	check_matlab
	
	forvalues test_price = $first_avgprice ($step_size) $last_avgprice {
		local test_price = round(`test_price',.1) //sometimes stata addes weird values like .000000000001 to the number in a loop
		global avgprice = "`test_price'"
		global avgprice = substr("${avgprice}",1,4) //still having issues sometimes... manually change the string...
		global RRnetwork = "_cfNORR"
		do "Theta_landregression_nobootstrap.do"
		* Now use that theta and avgprice pair to test the model estimated output and total shipments
		global best_theta = ${besttheta_AvgPrice}
		shell matlab -r "cd '${main_path}/structural';avgpricelist=[${avgprice}];thetalist=[${best_theta}];network=0;robustness_str='${robustness_str}';run Estimate_Theta_with_Total_Shipments.m;exit"
		* Wait for it to finish
		check_matlab
		* Note, this outputs the file: "output/theta_shipments${robustness_str}.csv" which
		* contains three columns, Theta;AvgPrice;Ratio of Total Shipments.
	}
	
	global firsttheta = 1.5
	global lasttheta = 6.5
	global iters = 100
	global first_avgprice = 37.9
	global last_avgprice = 39.9
	global step_size = .2

	shell matlab -r "cd '${main_path}/structural';fst=${firsttheta};lst=${lasttheta};iters=${iters};avgprice_list=[${first_avgprice}:${step_size}:${last_avgprice}];run Data_and_Theta_Prep_gridsearch.m;exit"
	* This waits for the last file to be created, should be VERY long
	check_matlab
	
	forvalues test_price = $first_avgprice ($step_size) $last_avgprice {
		local test_price = round(`test_price',.1) //sometimes stata addes weird values like .000000000001 to the number in a loop
		global avgprice = "`test_price'"
		global avgprice = substr("${avgprice}",1,4) //still having issues sometimes... manually change the string...
		global RRnetwork = "_cfNORR"
		do "Theta_landregression_nobootstrap.do"
		* Now use that theta and avgprice pair to test the model estimated output and total shipments
		global best_theta = ${besttheta_AvgPrice}
		shell matlab -r "cd '${main_path}/structural';avgpricelist=[${avgprice}];thetalist=[${best_theta}];network=0;robustness_str='${robustness_str}';run Estimate_Theta_with_Total_Shipments.m;exit"
		* Wait for it to finish
		check_matlab
		* Note, this outputs the file: "output/theta_shipments${robustness_str}.csv" which
		* contains three columns, Theta;AvgPrice;Ratio of Total Shipments.
	}
	
	import delimited using "output/theta_totalshipments${robustness_str}.csv", clear
	rename v3 sse
	twoway (scatter sse v2, mcolor(black)), ytitle("Squared Difference from Total Shipments") xtitle("Average Price of Transported Goods") graphregion(color(white))
	graph export "figures/fitting_total_shipments.png", as(png) replace
	qui sum sse
	* Keep theta and price that give the best fit to total_shipments. We have to
	* do this clunky adjustment because stata rounds when you grab the min
	keep if sse <= `r(min)' + 0.0000001
	count
	if `r(N)' > 1 {
		* If there are two tied for best fit, and thus fit the land regression,
		* they are near identical replicas (it is likely that a better fitting
		* theta and average price would be found in a smaller grid)
		keep if _n == 1
	}
	qui sum v1
	* Note no thetas are more than 4 decimals
	global globaltheta = round(`r(mean)',.0001)
	global globaltheta_baseline = ${globaltheta}
	qui sum v2
	global avgprice = round(`r(mean)',0.1)
	global avgprice_baseline = ${avgprice}

	* Calculate confidence interval with bootstrapping
	do "Theta_landregression.do"
}
else {

	global firsttheta = 1.5
	global lasttheta = 6.5
	global iters = 100
	* Since this branch of the if statement was created knowing the results, we hard code
	* the best option to reduce run time.
	global avgprice_baseline = 38.7
	global avgprice = ${avgprice_baseline}
	global avgprice_list = ${avgprice_baseline}
	global globaltheta_baseline = 3.05
	global globaltheta = ${globaltheta_baseline}

	* Launches Matlab and 
	shell matlab -r "cd '${main_path}/structural';avgprice_list=${avgprice_list};fst=${firsttheta};lst=${lasttheta};iters=${iters};run Data_and_Theta_Prep.m;exit"
	* This waits for the last file to be created, should be almost a full day
	check_matlab
	
	global RRnetwork = "_cfNORR"

	
	do "Theta_landregression.do"
	shell matlab -r "cd '${main_path}/structural';avgpricelist=[${avgprice}];thetalist=[${globaltheta}];network=0;robustness_str='${robustness_str}';run Estimate_Theta_with_Total_Shipments.m;exit"
	* Wait for it to finish
	check_matlab
	* Note, this outputs the file: "output/theta_shipments${robustness_str}.csv" which
	* contains three columns, Theta;AvgPrice;Ratio of Total Shipments.
	* The third row should be within 2 thousands of 1.
}

* Retrieve upper and lower bound of theta 95% CI, given price
	import excel using "tables/Theta_CI.xlsx", firstrow sheet("CI_AvgPrice${avgprice}") clear
	* Floor at 4 decimal places(for naming conventions) There are some inconsistencies
	* in rounding and machine precision between matlab and stata when we use loops.
	* They are irrelevant for the results as the changes are far to small. However,
	* for naming the files we need it to be precise. Thus, we just drop decimals
	* after 4 places.
	replace output_p2_5 = floor(output_p2_5*10000)/10000
	replace output_p97_5 = floor(output_p97_5*10000)/10000
	sum output_p2_5
	global lowtheta = `r(mean)'
	sum output_p97_5
	global hightheta = `r(mean)'

* Now we run many counterfactuals using the appropriate theta and average price
* Initialize a global for native born calculations
global nativeborn_str = ""
* First baseline, but run it for all years
global thetalist = "${globaltheta}"
* Restore the avg price to hold its original, baseline value
global avgpricelist = "${avgprice_baseline}"
* Various counterfactual networks, 0 -> no railroads ; 1850 -> 1850 RR network ; 1860 -> 1860 network ; 1870 -> 1870 railroads ;
* 1880 -> 1880 railroads ; 11 -> no railroads, include proposed canals ; 22 -> No RR, but cheaper wagons ;
* 33 -> No RR, canals and cheaper wagons ; 44 -> No RR, double waterway costs ; 55 -> No RR, canals and double waterways
* 66 -> double the cost of railroads
global networklist = "0"
* We want to run our main counterfactual (noRRs) for all years.
global yearlist = "1860 1870 1880 1890 1900"

shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
check_matlab

* this specifies whether to calculate CF MA for all decades. This is the only
* time we want to do that, so right after running this, we'll specify only_1890 = 1
global only_1890 = 0
global network = 0
* Install ereplace
cap ssc install ereplace


* Run Stata analysis
foreach year of numlist 1860(10)1900 {
	global year_of_interest = `year'
	do "Counterfactuals_analysis.do"
}

* Change capital assumptions--owned in NY
global yearlist = "1890"
global K_share_str = "NYC"
shell matlab -r "cd '${main_path}/structural';K_share_str='${K_share_str}';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals_ChangingCapital.m;exit"
check_matlab

global only_1890 = 1
global network = 0
* Install ereplace
cap ssc install ereplace

* Run Stata analysis
global year_of_interest = 1890
do "Counterfactuals_analysis_ChangingCapital.do"


* Capital Shares
use "output/model_implied_cfNORR_AvgPrice_38.7_data1880_theta3.05.dta", clear
keep fips year f_capital cf_capital_Uconst
egen tot_cap_f = total(f_capital), by(year)
egen tot_cap_cf = total(cf_capital_Uconst), by(year)
gen cap_share = f_capital / tot_cap_f 
gen cap_share_cf = cf_capital_Uconst / tot_cap_cf
keep fips year cap_share*
keep if year==1890
sort fips
merge 1:1 fips using "intermediate/Fips1890_cfNORR_data1880.dta", keep(match) nogen
drop year
preserve
	keep cap_share
	replace cap_share = 0 if cap_share == .
	export delimited "output/FACT_capital_shares_1890.csv", replace novar
restore
keep cap_share_cf
replace cap_share_cf = 0 if cap_share_cf == .
export delimited "output/CF_capital_shares_1890.csv", replace novar

use "input/wealth_shares.dta", clear
merge 1:1 fips using "intermediate/Fips1890_cfNORR_data1880.dta", keep(using match)
replace persprop_share = 0 if persprop_share==.
replace wealth_share = 0 if wealth_share==.
sort fips
preserve
	keep persprop_share
	replace persprop_share = 0 if persprop_share == .
	export delimited "output/PERS_capital_shares_1890.csv", replace novar
restore
keep wealth_share
replace wealth_share = 0 if wealth_share == .
export delimited "output/W_capital_shares_1890.csv", replace novar


* FACTUAL CAPITAL SHARES
global K_share_str = "FACT"
shell matlab -r "cd '${main_path}/structural';K_share_str='${K_share_str}';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals_ChangingCapital.m;exit"
check_matlab

global only_1890 = 1
global network = 0
* Run Stata analysis
global year_of_interest = 1890
do "Counterfactuals_analysis_ChangingCapital.do"

* COUNTERFACTUAL CAPITAL SHARES
global K_share_str = "CF"
shell matlab -r "cd '${main_path}/structural';K_share_str='${K_share_str}';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals_ChangingCapital.m;exit"
check_matlab

global only_1890 = 1
global network = 0
* Run Stata analysis
global year_of_interest = 1890
do "Counterfactuals_analysis_ChangingCapital.do"

* PERSONAL PROPERTY CAPITAL SHARES
global K_share_str = "PERS"
shell matlab -r "cd '${main_path}/structural';K_share_str='${K_share_str}';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals_ChangingCapital.m;exit"
check_matlab

global only_1890 = 1
global network = 0
* Run Stata analysis
global year_of_interest = 1890
do "Counterfactuals_analysis_ChangingCapital.do"

* TOTAL WEALTH CAPITAL SHARES
global K_share_str = "W"
shell matlab -r "cd '${main_path}/structural';K_share_str='${K_share_str}';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals_ChangingCapital.m;exit"
check_matlab

global only_1890 = 1
global network = 0
* Run Stata analysis
global year_of_interest = 1890
do "Counterfactuals_analysis_ChangingCapital.do"

global yearlist = "1890"
*HK Exercise
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals_HKexercise.m;exit"
* Wait for it to finish
check_matlab

global only_1890=1
global year_of_interest = 1890
do "Counterfactuals_analysis_HKexercise.do"

global only_1890 = 0
global networklist = "1850 1860"
global yearlist = "1860 1870 1880 1890 1900"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
check_matlab

foreach network in 1850 1860  {
	global network = `network'
	do "Counterfactuals_analysis.do"
}

* Do network counterfactuals
global only_1890 = 1
global year_of_interest=1890
global networklist = "11 66 1860 1870 1880"

global yearlist = "1890"

shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab

foreach network in 11 66 1870 1880 {
	global network = `network'
	do "Counterfactuals_analysis.do"
}

global robustness_str = "_data1880"

* AvgPrice = 20 | AvgPrice = 50
global networklist = "0"
global network = "0"
global avgpricelist = "20 50"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab

foreach avgprice_temp in 20 50 {
	global avgprice = `avgprice_temp'
	do "Counterfactuals_analysis.do"
}

* Theta = upper and lower bound on confidence intervals
global avgpricelist = "${avgprice_baseline}"
global avgprice = "${avgprice_baseline}"
global thetalist = "${lowtheta} ${hightheta} 8.22"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab

foreach theta_temp in $thetalist {
	global globaltheta = `theta_temp'
	do "Counterfactuals_analysis.do"
}

* No Agriculture Wedge
global globaltheta = ${globaltheta_baseline}
global thetalist = "${globaltheta}"
global robustness_strtemp = "_noagwedge"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_strtemp}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab

global robustness_str = "${robustness_strtemp}"
do "Counterfactuals_analysis.do"

* Native Born population
global robustness_str = "_data1880"
global nativeborn_str = "_nb"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab

do "Counterfactuals_analysis.do"


global nativeborn_str = "_nbparent"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab

do "Counterfactuals_analysis.do"

global nativeborn_str = "_true"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab

do "Counterfactuals_analysis.do"

global nativeborn_str = "_0"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab
do "Counterfactuals_analysis.do"


global nativeborn_str = "_10"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab
do "Counterfactuals_analysis.do"

global nativeborn_str = "_25"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab
do "Counterfactuals_analysis.do"

global nativeborn_str = "_75"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab
do "Counterfactuals_analysis.do"

global nativeborn_str = "_90"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab
do "Counterfactuals_analysis.do"

global nativeborn_str = "_100"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab
do "Counterfactuals_analysis.do"


global nativeborn_str = "_pop1860"
global networklist = "1860"
global robustness_str = "_data1880"
global yearlist = "1890"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"
* Wait for it to finish
check_matlab

global network = "1860"
do "Counterfactuals_analysis.do"

global nativeborn_str = ""
global networklist = "0"
global network = "0"
global robustness_str = "_KwedgeM"
global yearlist = "1890"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"

check_matlab
do "Counterfactuals_analysis.do"


global nativeborn_str = ""
global networklist = "0"
global robustness_str = "_KLwedgeM"
global yearlist = "1890"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"

check_matlab
do "Counterfactuals_analysis.do"

global nativeborn_str = ""
global networklist = "0"
global robustness_str = "_Kwedge0"
global yearlist = "1890"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"

check_matlab
do "Counterfactuals_analysis.do"


global nativeborn_str = ""
global networklist = "0"
global robustness_str = "_shrunk5"
global yearlist = "1890"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"

check_matlab
do "Counterfactuals_analysis.do"	

global nativeborn_str = ""
global networklist = "0"
global robustness_str = "_shrunk10"
global yearlist = "1890"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"

check_matlab
do "Counterfactuals_analysis.do"

global nativeborn_str = ""
global networklist = "0"
global robustness_str = "_shrunk25"
global yearlist = "1890"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"

check_matlab
do "Counterfactuals_analysis.do"

* Only Capital
global nativeborn_str = ""
global networklist = "0"
global robustness_str = "_shrunk5K"
global yearlist = "1890"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"

check_matlab
do "Counterfactuals_analysis.do"	

global nativeborn_str = ""
global networklist = "0"
global robustness_str = "_shrunk10K"
global yearlist = "1890"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"

check_matlab
do "Counterfactuals_analysis.do"

global nativeborn_str = ""
global networklist = "0"
global robustness_str = "_shrunk25K"
global yearlist = "1890"
shell matlab -r "cd '${main_path}/structural';yearlist=[${yearlist}];avgpricelist=[${avgpricelist}];thetalist=[${thetalist}];networklist=[${networklist}];robustness_str='${robustness_str}';nativeborn_str='${nativeborn_str}';run Run_Counterfactuals.m;exit"

check_matlab
do "Counterfactuals_analysis.do"
}
else {
global globaltheta_baseline = 3.05
global avgprice_baseline = 38.7
global avgprice = ${avgprice_baseline}
global globaltheta = ${globaltheta_baseline}
import excel using "tables/Theta_CI.xlsx", firstrow sheet("CI_AvgPrice${avgprice}") clear 
replace output_p2_5 = floor(output_p2_5*10000)/10000
replace output_p97_5 = floor(output_p97_5*10000)/10000
sum output_p2_5
global lowtheta = `r(mean)'
sum output_p97_5
global hightheta = `r(mean)'
}

* Make sure "model_MA_terms" file has all years for next steps
global only_1890=0
global yearlist = "1860 1870 1880 1890 1900"
global robustness_str = "_data1880"
global networklist = "0"
global nativeborn_str = ""

do "Counterfactuals_analysis.do"

********************************************************************************
***** GENERATE: POST-STRUCTURAL ************************************************
********************************************************************************

global avgprice = ${avgprice_baseline}
global globaltheta = ${globaltheta_baseline}
cd ../generate_poststructural
do "generate_poststructural.do"

********************************************************************************
***** REDUCED FORM RESULTS AND ROBUSTNESS **************************************
********************************************************************************

cd "${main_path}/reduced_form"
do "figures.do"
do "tables.do"

cd ..
